********************************************************************************
************************** DEMAND  EFFECT BOUNDS *******************************
********************************************************************************
*! version 1.0, Sep 26, 2016
*! author Johannes Haushofer
*! Demand Effect Bounds (de Quidt, Haushofer, Roth 2016)

capture program drop demandbounds
program demandbounds, eclass
	syntax varlist(min=1 max=2) [if] [in] [pweight fweight iweight/], [Demand(str)] [LEVel(real 10)] [CI] [COMmand(str)] [Robust]
	dis "`varlist'"
	tokenize `varlist'
	local y = "`1'"
	local t = "`2'"
	
	//capture confirm variable `2'
	if "`t'" == "" { // no treatment variable supplied; generate a dummy 
		tempvar t 
		gen `t' =1 
		label var `t' "Mean"
		local meanonly = 1
		dis "no t supplied. Meanonly: `meanonly'"
		}
		
//	quietly {
	* Choose regression command
	if "`command'" == "" {
			local command "reg"
		}
	
	* Set up "if"
		if "`if'" == "" {
			local if2 = "if"
		}
		else {
			local if2 = "`if' & "			
		}

	* Set level
		local clevel = `c(level)'
		if `level' >= 50 & `level' < 100 {
			local level = round(`level',0.01)
			set level `level'
		}
		else {
			set level `clevel'
			local level = `clevel'			
		}	
	
	tempvar d_ub d_lb 
	gen `d_lb' = (`t' == 0 & `demand' ==  1) | (`t'==1 & `demand' == -1)
	gen `d_ub' = (`t' == 0 & `demand' == -1) | (`t'==1 & `demand' ==  1)

	mat V = J(2,2,0)

	* Run estimation
	if "`meanonly'" ~= "1" {
		`command' `y' `t' `if2' `d_lb' `in', `robust'
		}
	else {
		`command' `y' `t' `if2' `d_lb' `in', `robust' nocons
		}
		ereturn list
		mat VV = e(V)
		mat V[1,1] = VV[1,1]
		local lb_b = _b[`t']
		local lb_n = `e(N)'
		local lb_se = sqrt(V[1,1])
		local lb_t = `lb_b'/`lb_se'
		local lb_df = `e(df_r)'
		local lb_p = (2 * ttail(`e(df_r)', abs(`lb_t')))

	if "`meanonly'" ~= "1" {
		`command' `y' `t' `if2' `d_ub' `in', `robust'
		}
	else {
		`command' `y' `t' `if2' `d_ub' `in', `robust' nocons
		}
		mat VV = e(V)
		mat V[2,2] = VV[1,1]
		local ub_b = _b[`t']
		local ub_n = `e(N)'
		local ub_se = sqrt(V[1,1])
		local ub_t = `ub_b'/`ub_se'
		local ub_df = `e(df_r)'
		local ub_p = (2 * ttail(`e(df_r)', abs(`ub_t')))
		
		mat b = J(1,2,0)
		mat b[1,1] = `lb_b'
		mat b[1,2] = `ub_b'
		
		matname b `t':lower `t':upper, columns(.) explicit
		matname V `t':lower `t':upper, explicit
	
	* Post results
	ereturn clear
	
	ereturn post b V, depname(`y') properties(b V) 

	ereturn scalar lb_n  = `lb_n'	
	ereturn scalar lb_b  = `lb_b'
	ereturn scalar lb_se = `lb_se'
	ereturn scalar lb_t  = `lb_t'
	ereturn scalar lb_df = `lb_df'
	ereturn scalar lb_p  = `lb_p'

	ereturn scalar ub_n  = `ub_n'
	ereturn scalar ub_b  = `ub_b'
	ereturn scalar ub_se = `ub_se'
	ereturn scalar ub_t  = `ub_t'
	ereturn scalar ub_df = `ub_df'
	ereturn scalar ub_p  = `ub_p'

	* Compute confidence interval for treatment effect
	if "`ci'" == "ci" {
//		if _rc == 0  {
			demandci
			ereturn scalar  lb_ci = cilower
			ereturn scalar  ub_ci = ciupper
			scalar drop ciupper cilower
		}
		else {
			ereturn local lb_ci "."
			ereturn local ub_ci "."
		}
//	}
	
//} // quietly 
	
	* Display results
	display _newline as text "Demand effect bounds (de Quidt, Haushofer, Roth, 2016)"
    display _newline as text "Number of lower bound obs." _skip(20) as text  " =   " as result `e(lb_n)'
    display			 as text "Number of upper bound obs." _skip(20) as text  " =   " as result `e(ub_n)'
	if "`ci'" == "ci" {
		display _newline as text "Confidence interval" _skip(20) as text  " =   [" as result `e(lb_ci)' as text "   " as result `e(ub_ci)' as text  "]"
    }
  	ereturn display, level(`level')
end


capture program drop demandci
program demandci, rclass
version 1.0
	local cs = invnormal(1-(100-`c(level)')/100) 
	local ce = invnormal(1-(100-`c(level)')/200)
	local qd = 10^15
	forvalues cc = `cs'(0.001)`ce' {
		local qdn = ((normal(`cc'+(_b[upper]-_b[lower])/max(_se[lower], _se[upper]))-normal(-`cc')) - (1-(100-`c(level)')/100))^2
		if `qdn' < `qd' {
			local qd = `qdn'
			local cnn = `cc'
		}
	}
	scalar cilower = _b[lower]-_se[lower]*`cnn'
	scalar ciupper = _b[upper]+_se[upper]*`cnn'
end

